//-----------------------------------------------------------------------------
//
// @file	Matrix.h
// @brief	Matrix  classes (code modified from DonHouse, Texas A&M Univ)
// @author	Vinod Melapudi
// @date	20Feb2010
//
//-----------------------------------------------------------------------------


#ifndef _MATRIX_H
#define _MATRIX_H

#include "Vector.h"


/* Matrix Descriptions and Operations */

class Matrix2x2;
class Matrix3x3;
class Matrix4x4;
class Matrix;

//
// The internal storage form of transforms.  The matrices in row-major
// order (ie  mat[row][column] )
//

class Matrix2x2 {
protected:
	Vector2 row[2];

public:
	Matrix2x2(float a11 = 0, float a12 = 0,
		float a21 = 0, float a22 = 0);

	Vector2& operator[](int i);
	const Vector2& operator[](int i) const;

	operator Matrix3x3();
	operator Matrix4x4();
	operator Matrix();

	void print(int w = 7, int p = 3) const; // print with width and precision

	void set(float a11 = 0, float a12 = 0,
		float a21 = 0, float a22 = 0);
	void identity();

	Matrix2x2 transpose() const;
	Matrix2x2 inv() const;

	friend Matrix2x2 operator+(const Matrix2x2& m1, const Matrix2x2& m2);
	friend Matrix2x2 operator-(const Matrix2x2& m1, const Matrix2x2& m2);
	friend Matrix2x2 operator*(const Matrix2x2& m1, const Matrix2x2& m2);
	friend Matrix2x2 operator*(float a, const Matrix2x2& m);
	Matrix2x2 operator*(float a) const;
	// mat times vector
	friend Vector2 operator*(const Matrix2x2& m, const Vector2& v);
	// vector times mat
	friend Vector2 operator*(const Vector2& v, const Matrix2x2& m);
	// outer product
	friend Matrix2x2 operator&(const Vector2& v1, const Vector2& v2);
};

class Matrix3x3 {
protected:
	Vector3 row[3];

public:
	Matrix3x3(float a11 = 0, float a12 = 0,
		float a21 = 0, float a22 = 0);
	Matrix3x3(float a11, float a12, float a13,
		float a21, float a22, float a23,
		float a31, float a32, float a33);

	Vector3& operator[](int i);
	const Vector3& operator[](int i) const;

	operator Matrix4x4();
	operator Matrix();

	void print(int w = 7, int p = 3) const;  // print with width and precision

	void set(float a11 = 0, float a12 = 0,
		float a21 = 0, float a22 = 0);
	void set(float a11, float a12, float a13,
		float a21, float a22, float a23,
		float a31, float a32, float a33);
	void identity();

	Matrix3x3 transpose() const;
	Matrix3x3 inv() const;

	friend Matrix3x3 operator+(const Matrix3x3& m1, const Matrix3x3& m2);
	friend Matrix3x3 operator-(const Matrix3x3& m1, const Matrix3x3& m2);
	friend Matrix3x3 operator*(const Matrix3x3& m1, const Matrix3x3& m2);
	friend Matrix3x3 operator*(float a, const Matrix3x3& m);
	Matrix3x3 operator*(float a) const;
	// mat times vector
	friend Vector3 operator*(const Matrix3x3& m, const Vector3& v);
	// vector times mat
	friend Vector3 operator*(const Vector3& v, const Matrix3x3& m);
	// outer product
	friend Matrix3x3 operator&(const Vector3& v1, const Vector3& v2);
};

class Matrix4x4 {
protected:
	Vector4 row[4];

public:
	Matrix4x4(float a11 = 0, float a12 = 0, float a13 = 0,
		float a21 = 0, float a22 = 0, float a23 = 0,
		float a31 = 0, float a32 = 0, float a33 = 0);
	Matrix4x4(float a11, float a12, float a13, float a14,
		float a21, float a22, float a23, float a24,
		float a31, float a32, float a33, float a34,
		float a41, float a42, float a43, float a44);

	operator Matrix();

	Vector4& operator[](int i);
	const Vector4& operator[](int i) const;

	void print(int w = 7, int p = 3) const;  // print with width and precision

	void set(float a11 = 0, float a12 = 0, float a13 = 0,
		float a21 = 0, float a22 = 0, float a23 = 0,
		float a31 = 0, float a32 = 0, float a33 = 0);
	void set(float a11, float a12, float a13, float a14,
		float a21, float a22, float a23, float a24,
		float a31, float a32, float a33, float a34,
		float a41, float a42, float a43, float a44);
	void identity();

	Matrix4x4 transpose() const;
	Matrix4x4 inv() const;

	friend Matrix4x4 operator+(const Matrix4x4& m1, const Matrix4x4& m2);
	friend Matrix4x4 operator-(const Matrix4x4& m1, const Matrix4x4& m2);
	friend Matrix4x4 operator*(const Matrix4x4& m1, const Matrix4x4& m2);
	friend Matrix4x4 operator*(float a, const Matrix4x4& m);
	Matrix4x4 operator*(float a) const;
	friend Matrix4x4 LU_Decompose(const Matrix4x4& M, int *indx);
	friend void LU_back_substitution(const Matrix4x4& M, int *indx, float col[]);
	// mat times vector
	friend Vector4 operator*(const Matrix4x4& m, const Vector4& v);
	// vector times mat
	friend Vector4 operator*(const Vector4& v, const Matrix4x4& m);
	// mat times vector
	friend Vector3 operator*(const Matrix4x4& m, const Vector3& v);
	// vector times mat
	friend Vector3 operator*(const Vector3& v, const Matrix4x4& m);
	// outer product
	friend Matrix4x4 operator&(const Vector4& v1, const Vector4& v2);
};

class Matrix {
protected:
	Vector *row;
	int Nrows, Ncols;

public:
	Matrix(int rows = 0, int cols = 0, const float *M = NULL);
	Matrix(const Matrix& M);
	Matrix(float a11, float a12,
		float a21, float a22);
	Matrix(float a11, float a12, float a13,
		float a21, float a22, float a23,
		float a31, float a32, float a33);
	Matrix(float a11, float a12, float a13, float a14,
		float a21, float a22, float a23, float a24,
		float a31, float a32, float a33, float a34,
		float a41, float a42, float a43, float a44);

	~Matrix();

	void setsize(int rows, int cols);
	int nrows() const;
	int ncols() const;

	operator Matrix2x2();
	operator Matrix3x3();
	operator Matrix4x4();

	Vector& operator[](int i);
	const Vector& operator[](int i) const;

	void print(int w = 7, int p = 3) const;  // print with width and precision

	void set(float *M);
	void set(float a11, float a12,
		float a21, float a22);
	void set(float a11, float a12, float a13,
		float a21, float a22, float a23,
		float a31, float a32, float a33);
	void set(float a11, float a12, float a13, float a14,
		float a21, float a22, float a23, float a24,
		float a31, float a32, float a33, float a34,
		float a41, float a42, float a43, float a44);
	void identity();

	Matrix transpose() const;
	Matrix inv() const;
	void svd(Matrix &U, Vector &W, Matrix &V) const;

	friend Matrix diag(const Vector &V);

	const Matrix& operator=(const Matrix& m2);
	friend Matrix operator+(const Matrix& m1, const Matrix& m2);
	friend Matrix operator-(const Matrix& m1, const Matrix& m2);
	friend Matrix operator*(const Matrix& m1, const Matrix& m2);
	friend Matrix operator*(float a, const Matrix& m);
	Matrix operator*(float a) const;
	friend Matrix LU_Decompose(const Matrix& M, int *indx);
	friend void LU_back_substitution(const Matrix& M, int *indx, float col[]);
	// mat times vector
	friend Vector operator*(const Matrix& m, const Vector& v); 
	// vector times mat
	friend Vector operator*(const Vector& v, const Matrix& m);
	// outer product
	friend Matrix operator&(const Vector& v1, const Vector& v2);
};

#endif
